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ABSTRACT 

We show how the 3DVAR data assimilation methodology can be used in the astrophys- 
ical context of a two-dimensional convection flow. We study the way this variational 
approach finds best estimates of the current state of the flow from a weighted average 
of model states and observations. We use numerical simulations to generate synthetic 
observations of a vertical two-dimensional slice of the outer part of the solar convec- 
tion zone for varying noise levels and implement 3DVAR when the covariance matrices 
are scalar. Our simulation results demonstrate the capability of 3DVAR to produce 
error estimates of system states between up to tree orders of magnitude below the 
original noise level present in the observations. This work exemplifies the importance 
of applying data assimilation techniques in simulations of the stratified convection. 



1 INTRODUCTION 



When using models to describe the temporal evolution of 
observed complex systems we are confronted with a number 
of challenges. An immediate difficulty in dealing with this 
question is that we generally do not know in all detail the 
current state of the system or the initial condition that is to 
be used. Lacking such information prevents us from keeping 
a model-based simulation in step with the behavior of the 
observed system. 

Data assimilation techniques offer means to address 
such challenges for complex systems by keeping a computer 
simulation (i.e. model) in synchronization with observations 
of the system it represents. It provides a general framework 
for simultaneously comparing, combining, and evaluating 
observations of physical systems and output from computer 
simulations. The methods used in data assimilation have 
been developed over several decades, primarily in meteo- 
rology and oceanography for the prediction of the future 
behavior. 

Data as similation is used daily in operational weather 
prediction dBengtsson et al.l Il98ll ). in climate forecast 
(|Palmer fc Hagedornll2006l ) and it was even used to correct 
the p ath of the A pollo spacecraft during the first moon land- 
ings (|Ciprall 19931 ). There is a large and growing body of lit- 
erature including s everal monographs (|Dalevlll993l ; iKalnavl 
120031 ; IWun sch 2006) a nd work discussing its theoretical foun- 
dations dLorend 198 ll ; lLorend[l98^ ; iLe Dimet fc Talagrandl 
1 19861 ; lGhillll989i r Astroph ysical data assimilation has re- 
cently been discussed by (jBrunl 12007). both in the con- 
text of space weather and in solar cycle prediction JPikpatil 
120071 : IChoudhuri et al.ll2007l ; iKitiashvili fc Kosovichevll2008l) 



as w ell as in the context of dynamo models ( Jouve et al.l 

[2011]). 



Here we focus on the three-dimensional variational 
(3DVAR) data assimilation technique, also known as sequen- 
tial approach |Dalcy 1993]), which produces updates of the 
current state of a model simulation at times when system 
observations are available. Propagation of model states be- 
tween times when the system is observed are free simulations 
of the model initiated at the latest state estimate. An exten- 
sion of 3DVAR to implicitly incorporate dynamical informa- 
tion is known as four-dimensional variational (4DVAR) data 
assimilation. 3DVAR dynamically evolves the mean state 
wheres 4DVAR also evolution other statistical properties of 
the model dynamics. 

State estimates produced by 3DVAR are optimal pro- 
vided that the model is linear and the uncertainties are 
Gaussian. In other words, 3DVAR states are Best Linear 
Unbiased Estimates, where best and optimal refe r to the 
lowest possible mean s quared error of the estimate (|Kalmanl 
Il960j ; iTalagrandl Il997h . 

For nonlinear models, error statistics may become non- 
Gaussian even when the initial distribution is Normal and 
3DVAR (or 4DVAR) estimates are not longer unbiased. In 
this case, data assimilation techniques are challenged by the 
fact that a ctual applications are typically based on nonlinear 
processes (jPires et al.ll 19961 ). Specifically, that the states ex- 
hibited by real systems under observation will diverge from 
those predicted by a model simulation is clear and this is 
principally owing to two causes (Palmer & Hagedorn 2006): 
observational error and sensitivity to initial conditions. The 
first of these is a result of what may be called noise. Since 
its statistical character may not be known, we may need to 
make some assumptions about its properties. The second 
source of error occurs on many complex systems and is re- 
ferred to as chaotic behavior. This has been known for some 
time, but only in recent decades serious progress in its un- 
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derstanding has been possible. Sensitivity of the model to 
initial conditi ons limits ho w far into the future predictions 
can be made l|Lorenj|l993l ). 

Among other challenges present for data assimilation of 
nonlinear model states like the mismatch between spatial lo- 
cations of observations and grid positions of the model, and 
the unpaired model variables to observed quantities make 
the study towards more effective data assimilation tech- 
niques an important and demanding area of research. De- 
spite these challenges and open questions, 3DVAR is widely 
used in the oceanographic and meteorological communities, 
and would make a good candidate method to be explored in 
the context of astrophysical flows. 



2 STRATIFIED CONVECTION MODEL 

We are motivated to use data assimilation techniques in the 
context of stratified convection as a path to obtain pre- 
dictions of solar subsurface weather events, i.e. the flow 
structure beneath the surface. Anticipating the possibility 
of violent events on the solar surface such as coronal mass 
ejections that affect the space weather and th e dynamics of 
the E arth's magnetosphere is important, (see lllonidis et al] 
2011). The idea is to use a model of solar subsurface convec- 
tion, ultimately involving the magnetic fields that give rise 
to surface activity such as coronal mass ejections, although 
curre nt attempts in that di rection are still at a preliminary 
stage (|Warnecke et alj|201ll ). However, once such models are 
able to reproduce sufficient details of solar activity, it will be 
important to synchronize the model with daily observations 
to be able to use it for predictions. 

As a proof of concept, we design a data assimilation ex- 
periment to test the implementation of 3DVAR for the Pen- 
cil Code, a public domain code of high-order (sixth order in 
space and third order in time) for solvin g the hydrodynamic 
equations l|Brandenburg fc DobleJ l2002n We consider here 
a simple two-dimensional convection model representing the 
turbulent flows of stars with outer convection zones. In our 
experiment, synthetic observations are generated by adding 
noise to the output from our model. These observations are 
then processed by 3DVAR to produce an analysis. An anal- 
ysis is an estimation of t he unknown state of a system in 
terms of model variables l|Lorendll986l ; lTalagrandlll997h . 

Our implementation is general and can in the future be 
used for other problems that can be addressed with the Pen- 
cil Code. In this work we assume the model to be ideal and 
reproduces the same features present in the observations. In 
real world applications, the models are far from ideal, and 
imperfections and uncertainties related to the model are al- 
ways present. Ideally, we would like to be able to account for 
some portion of those unknowns by using data assimilation 
techniques. 

We use the sample 2d/conv-slab-MLT of the Pencil 
Code (revision rl4696 and later). This sample simulates a 
vertical two-dimensional slice of the outer part of a stellar 
convection zone. In particular, we use it to simulate con- 
vection at low resolution, 64 x 64, at a Rayleigh number of 
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8 x 10 5 l|Dintrans et al.ll2005l ). and a Reynolds number of ap- 
pr oximately 30. The basic se tup is similar to that described 
in dBrandenburg et al.|l2005l) and other models before them 
l|Hurlburt et al.lll986l ; iBrandenburg et al]|l996l ). consisting 
of a convectively unstable layer sandwiched between two sta- 
ble layers. 

The simulated vertical two-dimensional slice of the 
outer part of the solar convection zone has a mean field ve- 
locity of M rms = 0.08, the wavenumber of the energy-carrying 
eddies is k{ = 2nd for a depth of the unstable layer d = 1, 
therefore the correlation time r cor = (u rms fcf) _1 is approxi- 
mately 2. 

Starting from an initial velocity field of perturbations 
with an amplitude of 3 x 10 -4 times the average sound speed, 
convective motion is generated without having to introduce 
any stochastic elements. This model is chosen to illustrate 
3DVAR in an astrophysical context for its sufficiently com- 
plex behavior without having any stochastic elements. 



3 DATA ASSIMILATION SETUP 

The 3DVAR scheme was developed in the meteorological 
community to improve model-based weather prediction in 
spite of observational and modeling uncertain ties. It was for - 
mulated in a unified Bayesian framework by (|Lorendll986l ). 
3DVAR produces updates of the current state of the system 
at times when observations are available, which in turn can 
be used as a new initial condition to be propagated forward 
to the time when the next observation is available. We can 
use 3DVAR as a black box along with a low resolution simu- 
lation to assimilate many data points at low computational 
cost on a laptop computer. For example, a typical model 
run for a 64 x 64 two-dimensional convection field over a 
time interval of 300 time units takes about 15 minutes on a 
laptop computer. 

3DVAR minimizes the sum of the squared differences 
between both the model background state and the observa- 
tions to find a solution that is a compromise between these 
two estimates of the true state. 

It is important to realize that in real problems, the true 
state is available only through noisy observations of the sys- 
tem. Therefore, it is impossible to tell how close our model 
output is to current a nd future true states of the system. In 
our twin- experiment (|Bengtsson et all 11981) we select two 
different initial fields to run the Pencil Code simulation. 
One of these initial fields represents the unknown true initial 
state of the system. The other initial field represents what 
might be, in practice, a good approximation or guess, of the 
initial state of the system. 

As mentioned before, we are set at the ideal case where 
there is no model uncertainty and the only source of uncer- 
tainty is in the observations. In this way, we can assess how 
far/close the model state is to the true state of the system. 
The key is to generate a known true state against which 
the estimated state obtained via data assimilation can be 
verified. 

The experiment is setup as follows: the initial field cho- 
sen to represent the true initial state of the system, initializes 
a model run considered to be the original state of the sys- 
tem to be used as reference trajectory or control. The other 
initial field is used to initialize two different runs: one a free 
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model run and other that will become the assimilated tra- 
jectory or analysis. The analysis is a collection of segments 
of model trajectories initialized at the 3DVAR corrections 
made at all times where the observations are available. The 
model (analysis) states at the time when an observation is 
available are known as the background states. 

Comparison of the free and the control run gives us a 
measure of the sensitivity to initial conditions of our model, 
i.e. it shows how similar initial conditions diverge in time. 

Similarly, comparing the free run and the analysis rep- 
resents the effect of the data assimilation procedure over a 
trajectory starting at the same initial field. If the assimila- 
tion of the second set of initial conditions is effective it will 
bring the analysis "closer" to the control run, and further 
away from the free run. 

We generate synthetic observations by adding indepen- 
dent and identically distributed noise to the horizontal ve- 
locity field in all grid-points to the control run. These syn- 
thetic data are considered to be our experimental observa- 
tions which in turn will be used to update the analysis. As 
is explained in the next section, 3DVAR requires both a 
model and observation state to find update the analysis at 
each given time in the assimilation window. 



4 3DVAR AND THE WEIGHT FACTOR 

The 3DVAR technique finds a model state X that agrees 
with the current state of the system and the information 
available in the observations and the model. Specifically, we 
minimize the weighted average of the residues for both ob- 
servations Yo and model states X b at time t to find an opti- 
mal solution. This is expressed in the 3DVAR cost function 
!|Lorendll986h 



J(X) =- [X - X b ] T B 1 [X - X b ] 

+±[Y -H(X)] T R- 1 [Y -H(X)]. 



(1) 



where Xb is the model state-traditionally called the back- 
ground state, Yo is the observed state, B is the background 
covariance matrix, R is the observational covariance matrix, 
and H is the observation operator. 

The analysis is X a , minimizes J(X) and corresponds 
to the best estimate of the current state of the system. At 
that time, the model is integrated using as initial condition 
the analysis up to the next time an observation is available. 

Synthetic observations are denoted by Yo, these obser- 
vations contain noise of amplitude, or, proportional to the 
maximum amplitude of the full 2D vertical velocity field. 
For example, a noise level of 1% corresponds to <tr = 0.01 
for a normalized field or or = 6 x 10~ 3 for an unnormalized 
field. 

The selection and construction of the observational and 
background covariance m atrices (R y B) i s of great inter- 
est in data assimilation (|Bannisterll2008al lbl). In our case, 
the observational noise is not correlated in space, and we 
neglect spatial correlations between model states, making 
the off-diagonal components of the B and R matrices van- 
ish. By doing this, we can set these matrices to be scalars, 



Ri 



SijQ-R and Bi 



Sij<r B . Without these spatial cor- 



less smooth over the two-dimensional domain. In more so- 
phisticated formulations of equation {TJ, the form of B can 
also include physical constraints to processes not resolved in 
the model. Literature in this ar ea is extensive, especially in 
oceanography and meteorology (|Dobricic fc Pina rdi 2008). 

In turn, we set the observation operator to be Hij = Sij. 
This means that we assume observations cover all grid- 
points in the model domain, i.e. the observables and model 
variables belong to the same space. In other words, system 
and model are the same. In more realistic implementations 
of 3DVAR, H is typically a computer algorithm that can- 
not be expressed explicitly as a matrix due to its nonlinear 
nature (|Dobricic fc Pinardil l2008). For example the transfor- 
mation between observables and variable might require the 
modeling of dynamical processes or making averages. 

After these assumptions over the matrices in Q, we can 
translate the cost function to: 



J(X) =w(X- X b f + {X- Yo) 



(2) 



where w is the ratio of the scalar variances corresponding to 
the observed and background states 



= (or/ob) 



(3) 



In this setting, to find the state v ector X a that min imizes 
we use Powell minimization (Pre ss et ai1ll992T ). 
The coefficient w in equation ([2]) behaves as a weight 

in the optimization process and will be referred to as the 

weight factor. 

Solving VJ(X) = yields 



w(X - X b ) + (X- Y ) 



0. 



(4) 



and the optimal state of the model that represents the sys- 
tems given the current observations, or analysis X a = X, is 
given by: 



X a = 



-Y + 



-X h 



(5) 



relations 3DVAR generates an analysis that is, in general, 



Equation ((5J unveils how the contribution of the background 
states X b and observations Yo affects the analysis X a in 
terms of the weight factor w. 

Figure [T] clearly shows the result of this process at times 
where observations (grey o) are available, 3DVAR performs 
a correction given by equation ((2| to a value referred to as 
analysis (black +), and from which a segment of background 
states (dotted curve) is initialized and run up to the next as- 
similation time. The analysis is the union of the background 
states for times different from the assimilation times, and 
the corrected values obtained at assimilation times. 

Given a model and a fixed set of observations, equation 
([5]) help us understand the effects of the weight factor in 
the resulting analysis X a , as it is detailed in the following 
paragraphs. 

For w < 1 or or < or, is interpreted as the case where 
we assumed that the observational uncertainty is smaller 
than the model uncertainty, weight is given to the obser- 
vations since small w will allow the distance (X a — -X";,) 2 
to grow without making large contributions to the cost 
function. In contrast, having a W > 1 will favor model 
states reflecting our assumption that there is more uncer- 
tainty related to the observations than that of the model, or 
or > OR. 

The next section presents and describes the results of 
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Figure 1. Temporal evolution of the horizontal velocity at a cer- 
tain midpoint of a 2D convection field for t £ [0, 500] . Control, 
free run, and analysis correspond to the grey, black, and dotted 
curves. A grey 'o' represent observations, and a '+' corrections 
(X a ) made by 3DVAR, both at assimilation time. 



our numerical experiments. We measure the "quality" of the 
obtained analyses using 3DVAR by varying values of the 
weight factor w as well as a couple of sets of observations 
with different noise levels. 

The correlation time was found to be approximately 2 
from simulation parameters (see Section [2}. This determine 
the relevant time scales of the convection features of our 
simulation and then we choose the most appropriate assim- 
ilation time. Consequently, our choice of assimilation time 
is large enough to let the oscillations propagate over the 
two-dimensional field but still small enough to be able to 
capture the smaller scale dynamics. In each 3DVAR exper- 
iment, data assimilation corrections are made each 10 time 
units (sound-travel-time units) and found no fundamental 
difference when using assimilation times between 5 and 20 
time units. 




300 




300 



Figure 2. Data assimilation run over observations marked with 
grey 'o' with 1% noise for w = (upper panel) and w = 0.5 (lower 
panel). Black '+' marked the 3DVAR correction at assimilation 
time. Grey, dotted, and dashed curves correspond to the control, 
analysis, and free trajectories, respectively. Both panels include 
zoom-ins for t e [170, 210]. 



5 RESULTS 

We generate analyses using 3DVAR for each of the weight 
factors w = {0,0.1,0.5,1,10}. For a fixed value of w we 
generate an analysis which is the result of assimilating one 
set of observations with either 1% and 2% noise level for the 
same set of two initial field conditions. 

The resulting horizontal velocity, u x , at the midpoint of 
the upper right quadrant of the two-dimensional domain, is 
plotted in Figure [2] for w — in the upper panel and w = 0.5 
in the lower panel, and in Figure [3] for w = 1 in the upper 
panel and w = 10 in the lower panel. 

Note that grey and dashed lines are the same in all 
panels since they represent the reference states of the system 
(control) and the corresponding free run of the model. 

Observations are plotted with grey 'o' symbols. Black 
'+' marks are used for the corrections calculated at assim- 
ilation time by minimizing expression Between assimi- 
lations, the analysis is composed by segments of model tra- 
jectories (dotted segments) initialized at the corrected state 
'+' as seen in detail at all insets in Figs. [5] and [3] From 
equation ([5]) and these insets, we clearly illustrate the am- 
plitude of the correction made in each case for each value 



of the weight factor w. This amplitude is measured by the 
gap between the '+' and the last dot of the previous dotted 
segment (background states). 

For the trivial case of setting w = 0, the first term in 
((2| is neglected and the best estimate of the current state of 
the system is given by X a ~ Yo from in (J5J, and as seen in 
the upper panel of Figure[2] The correction ('+') is "pulled" 
from the background state (dotted curve) to the observa- 
tion ('o') at assimilation time. The background states are 
just segments of transient trajectories initialized at the ob- 
servations. 

For any other value of w the starting point of the cor- 
rection is in between the observation ('o') and the end of 
the previous background states segment (dotted lines); see 
insets in Figs. [2] and [3] 

For < w < 1, the optimal value of the cost-function 
© is by a factor 1/w closer to the observations, Y , than 
to the background state Xt; see Eq. ([5}. The lower panel 
of Fig. [2] shows results for w — 0.5 where the corrections 
('+') fall closer to the observations than to the end of the 
last segment of background states. This figure shows that 
the analysis follows the control trajectory closely (solid grey 
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Figure 3. Data assimilation run over observations marked with 
grey 'o' with 1% noise for w = 1 (upper panel) and w = 10 
(lower panel). Black '+' marked the 3DVAR correction at as- 
similation time. Grey, dotted, and dashed curves correspond to 
control, analysis, and free trajectories, respectively. Both panels 
include zoom-ins for t £ [170,210]. 



line). Note that even if the corrections are large, e.g. at t = 
50 or t — 90, the analysis quickly relaxes to the control state. 

In the case w = 1, equal weights are given to model 
states and observations. The optimal value of X a is the av- 
erage of Yo and Xb, from ((SJ). No preference is given to any 
estimate and the midpoint is the optimal choice for X a as 
seen in the upper panel of Figure [3] 

When w > 1, the optimization of equation ([2]) will favor 
model states rather than observations as follows from equa- 
tion (0). The analysis at assimilation times is by a factor w 
closer to the model states than to the observations. In the 
lower panel of Fig. [3] we plot the resulting trajectories for 
this case. 

For w = 10, we observe from this plot and from results 
at other locations of the two-dimensional domain that for 
w > 1 the estimates of the original state of the system are 
biased toward the background states. 

When the model is linear and ideal, the analysis pro- 
duced by 3DVAR will be an unbiased estimate of the true 
state of the system. There, information about the system is 
only contained in the experimental observations. When the 
model is not ideal nor linear, information about the known 



unknown processes important to represent the system of in- 
terest but not part of the model, could only be included via 
the w. In more general terms, through the covariance ma- 
trices R and B which informed the 3DVAR procedure about 
model and observational known uncertainties. 

In this simplified experiment, we can see the great im- 
portance of the weight factor w in 3DVAR. It points out how 
crucial the construction of the covariance matrices B and R 
are for optimal results when using variational approaches of 
data assimilation. One of the motivations for our choice of 
R and B to be scalars is to set a baseline from which we can 
illustrate, in a simplified way, the inner workings of 3DVAR. 
It can be hard to see how the different components interact 
to create a result when more sophisticated choices of R and 
B are used. 

We would like to note that the zoom-in is chosen to be 
in the interval t £ [170, 190], in all panels of Figs. [2] and |3] 
as an example of an interval where 3DVAR does not per- 
form very well, and systematically pulls the analysis away 
from the control trajectory - considered here as the original 
system trajectory. The reason for this behavior calls for fur- 
ther studies but it is worth noticing that around t = 200 the 
performance of the assimilation returns to its previous level. 
The case w > 10, presented in the lower panel of figure 3, 
actually performs better during this interval. 

Consistently, we observe that the analysis is on average 
closer to the control trajectory than the observations for 
all values considered for w < 1. As noted, exceptions are 
observed for larger values of the weight factor and during 
the interval shown in the insets of Figure 2 and 3. 

Tables [1] and [2] present several measures of variation of 
the output from the simulations of our twin experiment. For 
each simulation, we calculate the variance of the distances 
between the control and the observations (Table [1]) or the 
control and the analysis (Table [5} over the data assimilation 
window. 

Specifically, Table [1] shows the variance of the distance 
between the control trajectory and the free trajectory (first 
row), and the noisy observations for 1% (second row) and 2% 
(third row) noise levels at the midpoint of the field (second 
column) denoted by (((8u®) 2 }} , this column measures the 
variability of the local behavior. In addition, the averaged 
variance over the whole vertical two-dimensional field for the 
first and the second half of the assimilation window (third 
and forth column), denoted by (((5 u x ) 2 }}t where the sub- 
script T = 1,2 refers to averages over the first or second 
half of the assimilation window, t £ [1, 150] or t £ [151, 300] 
respectively. 

Note that from the values in Table [TJ the free run is one 
or two orders of magnitude further away from the original 
state of the system (control) than the 1% and 2% noisy 
observations. Large difference between the second and third 
columns reflect how free run is diverging from the control 
run over two different time intervals. 

Table [1] is the baseline from which we measure the per- 
formance of 3DVAR when estimating the system state from 
noisy observations or the free trajectory-or 0% noise level-. 
Note that the variance of the noise is of the order of 10~ 6 . 
When assessing the performance of 3DVAR for different val- 
ues of w and noise levels, we look for variability measures 
between the control and the analysis lower than the levels 
set by the free run and the noisy observations. 
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Noise Level ((fag)*) ({(Sug) 2 ))! {{(5ug) 2 )) 2 

free 1100 530 3000 

1% 27 40 36 

2% 110 140 140 



Table 1. Measures of variability (in 10 -6 ) for the distance be- 
tween the control and observations with noise levels 1%, 2%, and 
the free trajectory. See text for notation description. 

All variance measures for the free trajectory are larger 
than the corresponding variances values for all w and noise 
levels. This means that performing 3DVAR data assimila- 
tion is more effective at estimating the original state of the 
system than just using as estimate the trajectory initialized 
with a very close initial field. The variance of the distance 
between the two initial velocity fields is w 2 x 10~ 7 . 

Table [2] shows the average variance of the distance be- 
tween control run and the analysis at the midpoint of the 
field (second column) denoted by (((Su^) 2 )) , and the aver- 
aged variance over the whole vertical two-dimensional field 
for the first and the second half of the assimilation win- 
dow (third and forth column), denoted by (((<5 a u x ) 2 ))t, for 
T = 1,2, in analogy of the notation used in Table [T] 

From the values presented in both tables, we observe 
that for both noise levels studied and all cases with w < 10, 
the variances of the distance between the control and the 
analysis are smaller-on average-than the corresponding ob- 
servations with the same noise level. This shows that on 
both local and global scales, over the two-dimensional do- 
main, 3DVAR is effectively reducing noise and accounting 
for sensitivity to initial conditions i.e. it is estimating a value 
for the horizontal velocity closer to the original state of the 
system (control) than to the trajectory generated using a 
guess of the initial state of the system (free run). 

In contrast, for w = 10 and both noise levels, the free 
and the analysis are-in average-consistently farther away 
from the control than the original noisy observations (corre- 
sponding values in Table fj} both locally and globally. This 
is also observed, for example at the specific location of the 
two-dimensional domain shown in figure [2] and [3] 

We note that during the time interval t G [170, 210] (see 
insets of Figs. [2] and [3]) this simplified version of 3DVAR 
under-performs compared to other times. More sophisti- 
cated choices for R and B, that inform more realistically 
the cost-function ((2j about spatial correlations over the field 
might generate an improved and consistent performance. 
When comparing in this way the variance of the distance 
between the control and the free run we write (((S f u :c ) 2 }}t, 
for T = 1, 2. 

Furthermore, in Figure]?] we plot in a semi- logarithmic 
scale the variance of the point-by-point distances of the two- 
dimensional vertical slices between the control and analysis 
(grey '•'), the control and the observations (grey 'o'), and 
the control and the free run (black '•') for all t G [1,300]. 
The black '+' mark the variance of distances between the 
control and analysis fields at assimilation times, i.e. when 
the corrections from ([5]) are made. 

From top to bottom, plots correspond to w — 0.1, 1, 10. 
Black dotted and grey curves are the same for all panel 
in the figure. The noise level is constant corresponding to 
1%. 



w 






«(<5u A ) 2 » 2 


Noise level 1% 





0.61 


0.34 


3.80 


0.1 


0.67 


0.35 


4.40 


1 


3.30 


4.10 


130.00 


10 


300.00 


500.00 


1600.00 


Noise level 2% 





2.40 


1.30 


19.00 


0.1 


2.70 


1.40 


22.00 


1 


33.00 


9.80 


840.00 


10 


950.00 


430.00 


2800.00 



Table 2. Measures of variability (in 10 — 6 ) for the distance be- 
tween the control and the analysis for several values of the w of the 
a;— component of the velocity, u x . The italic values are larger then 
the corresponding noise level. See text for notation descrition. 



It is important to note, that values plotted in Figure [4] 
are not the running variance of the distance between con- 
trol trajectories and the other relevant trajectories but the 
variance of the distances between the 2D fields at each time 
step. 

We observe in Fig. [4] that for all values of w and for t G 
[0, 10] the analysis and free runs are good estimators of the 
control trajectory. In addition, for all w and t G [20, 100], the 
free run and the analysis variances with respect to control, 
are below the noise level, with the grey curve below the black 
curve, showing that for 1% noise level 3DVAR produces a 
better estimate than the free run. Only for w = 0.1 is this 
the case all values of t G [30, 300]. Note the increase in the 
error of the estimate from t G [160, 190]. This is part of the 
interval shown in the insets for Figs. [2] and [3} . 

Otherwise, the the error of the estimate exceeds the 
noise-level for t G [170, 300] for w = 1, and for t G [190, 230] 
for w = 10. 

Figure [4] also shows how far apart-in average-the dis- 
tance between the control and the free run grows as time 
increases. The increase in the amplitude of the variance (the 
black curve) is up to six or seven orders of magnitude over 
the assimilation window with respect to the same distance 
at small times. As noted before, this exhibits the mode! sen- 
sitivity to initial conditions and it is also a feature observed 
for all values of the weight factor and noise levels. 

In the top panel of the Figure [4] we again see evidence 
of how the 3DVAR correction (black '+') lands closer to 
the observations and how the background states start to 
move towards to the control trajectory for lower values of 
the weight factor, w. The effect of w < 1 in this case, is to 
"pull" the background states back closer to the observations 
as the grey dotted curve is below the corrections ('+')• 

We recall that the weight factor is w = (or/ob) 2 , and 
it is interpreted here as a measure of the relative confidence 
given to either the observations and the model. High values 
of w reflect higher trust in the observations than in model 
representation of the system, and the opposite is valid for 
low values of w. Note that ctb relates to initial condition 
sensitivities than with model deficiencies in our simple set- 
ting. As commented earlier, the estimate is sensitive to the 
choices of ere and a more sophisticated choice can include 
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Figure 4. Semi-logarithmic plot of the variance of the point- 
by-point distances of the two-dimensional vertical slices between 
the control and the analysis in grey '•' (((8 A \i x ) 2 )), the control 
and the observations in grey 'o' (((5° and the control and 
the free run in black '•' (((5 F u x ) 2 )) for all t G [1,300]. The 
black '+' corresponds to {(8 A u x )' 2 ) at assimilation times, after 
the correction is made. All quantities are scaled by 10 6 . 

addition components that might help alleviate some of the 
model deficiencies. 

In conclusion from the results of our experiment, in the 
particular case where the covariance matrices are assumed 
to be scalar and for small values of the weight factor w, 
the effect of trusting the observations more than the model 
states (ur < <tb), provides a closer estimate of the original 



state of the system than just generating a trajectory close 
enough in initial conditions. This means that 3DVAR is suc- 
cessful at finding an optimal state estimator in the limit of 
small observational noise, model uncertainty related only to 
sensitivity to initial conditions, and scalar covariance matri- 
ces, R and B. Our results strongly point out the importance 
of choosing more sophisticated covariance matrices (R and 
B) to better inform the assimilation procedure about the 
known uncertainty sources in a problem of interest. 



6 DISCUSSION 

We have presented an idealized case where model and system 
are the same: A computer simulation is used both to gener- 
ate synthetic observations and as the model required for the 
data assimilation procedure. In this way, we can assess how 
far/close the model state estimates are to the true state of 
the system. The key is to have access to the true states that 
we can use to verify and evaluate estimates obtained using 
data assimilation. 

We used a simplified formulation of the 3DVAR data 
assimilation technique in terms of the weight factor: w = 
{<jr/ <jb) 2 , that defines the contribution of the model 
states -that contain propagated information from previous 
observations- and the current observation to make a state 
estimate. 

This formulation of 3DVAR used here is achieved by 
reducing the covariance matrices, R and B, to scalars; and 
the observation operator, H, to the identity. This selection 
corresponds to neglecting all spatial correlations between 
model states over the two-dimensional domain in addition to 
one-to-one correspondence between system observables and 
model variables. In this way, we clearly separate the con- 
tribution of observations and model states to the estimated 
state, as seen in equation ((2| and ([5]). 

It is less direct to see how the different components in- 
teract to create an estimate of the original state of the sys- 
tem when more sophisticated choices of covariance matrices, 
that represent uncertainties and spatial correlations. In that 
case we would have to think about the optimal combina- 
tion in analogy to equation ([5]), in terms of a generalization 
of the weight factor w as a weight matrix, W = RB _1 . In 
this analogy, model states and observations will be projected 
by the matrices W[W + I] -1 and [W + I] -1 , respectively, to 
contribute to the analysis X a . Here, I is the identity matrix. 

In general, we can say that to understand the 3DVAR 
algorithm it is important to look at the weight factor, par- 
ticularly in the limit where W is assumed to be the scalar 
w. 

Consistently we observe that the error between the state 
estimate and the original state of the system is below the 
noise level when more weight is given to the observations 
then to the model state. When the contribution from the 
model state is larger than the contribution from the obser- 
vations we note that the error eventually becomes larger 
than the noise level, see case w = 10. 

We note in Figure [2] [3] and H that 3DVAR under- 
performs for t G [160, 230] both locally and globally for all 
values of w. Further study of the simulation is needed to 
account for this atypical behavior. 

Minute differences in initial conditions generate differ- 
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ent time evolution for the different runs, as is expected for 
nonlinear systems. This is illustrated by the black curves 
in Figure [4] that present how the variance of the distance 
between the two initial conditions grows over the time in- 
terval. It can also be seen in Figure Q] where the grey and 
black curves, that started with close initial condition, are 
very different at latter times. 

On the other hand, large correction made by 3DVAR, 
for example the black '+' at time 50 and 90 in Figure [21 
does not appear to have a strong effect. The model run that 
starts at these far away estimates converges almost instantly 
back to the control at the same time. 

These features, that is reminiscent of chaotic behavior, 
can be understood in terms of attracting sets where small 
changes in initial condition generate a different time evolu- 
tion on the attracting set. A large correction probably takes 
us outside the attracting set and the solution rapidly falls 
back when the model is integrated forward. 

Another interpretation is that the large correction takes 
us to states that are not consistent with conservation laws 
and other physical constraints. The system would then 
rapidly be forced back onto a more physical state. This 
dual picture, using both physical and mathematical intu- 
ition somewhat clarifies this contradictory behavior. 

The 3DVAR methodology is optim al for linear m odels 
and Gaussian distributed uncertainties (|Lorendll986l ). Very 
few models in astrophysics have these properties. The va- 
lidity of variational methods outside the linear or weakly 
nonlinear case is unclear but if the assimilation is frequent 
enough the behavior might be closer to linear. Higher assim- 
ilation frequency during the interval t £ [170, 230] might for 
example give a better result. 

The indication of chaotic properties and attracting 
state space sets invites the use of other data assimilation 
methods that explic i tly take these pr operties into account 
jjudd fc Smithll200ll ; [judd et alj200Sh . They might be more 
applicable to non-linear astrophysical processes. 

In this work we have used an ideal model, that is, the 
system and the model is one and the same. There will al- 
ways be some limitation to modeling of real system and this 
might prove problematic. In general it is, in our mind, more 
important for models used for prediction of real systems to 
be reasonably realistic rather than just being close to the 
observations. 

For nonlinear astrophysical systems that operate on 
timescales from seconds to years, data assimilation will be 
of fundamental importance when quantitative agreement be- 
tween model and observations is to be assessed. The ultimate 
verification that a model is correct is its ability to make re- 
liable prediction and for this, data assimilation is necessary. 
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